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We present a finite element algorithm that computes eigenvalues and eigenfunctions of the Laplace 
operator for two-dimensional problems with homogeneous Neumann or Dirichlet boundary condi- 
tions or combinations of either for different parts of the boundary. In order to solve the generalized 
eigenvalue problem, we use an inverse power plus Gauss-Seidel algorithm. For Neumann boundary 
conditions the method is much more efficient than the equivalent finite difference algorithm. We 
have cheked the algorithm comparing the cumulative level density of the espectrum obtained numer- 
ically, with the theoretical prediction given by the Weyl formula. A systematic deviation was found. 
This deviation is due to the discretisation and not to the algorithm. As an application we calculate 
the statistical properties of the eigenvalues of the acoustic Bunimovich stadium and compare them 
with the theoretical results given by random matrix theory. 



Presentamos un algoritmo de elementos finitos que calcula eigenvalores y eigenfunciones del operador 
de Laplace para problemas en dos dimensiones con condiciones a la frontera de Neumann o Dirichlet o 
combinaciones de ambas en distintas partes de la frontera. Para resolver el problema de eigenvalores 
generalizado, usamos un algoritmo de potencias inverso mas otro de Gauss-Seidel. Para condiciones a 
la frontera de Neumann, el metodo es mucho mas eficiente que el algoritmo equivalente de diferencias 
finitas. Hemos probado el algoritmo comparando la densidad acumulada de niveles del espectro 
obtenido numericamente, con la prediction teorica dada por la formula de Weyl. Se encontro una 
desviacion sistematica. Esta desviacion es debida a la discretization y no al algoritmo. Como una 
aplicacion, calculamos las propiedades estadfsticas de los eigenvalores del estadio de Bunimovich 
acustico y las comparamos con los resultados teoricos dados por la teon'a de matrices aleatorias. 

Subject Classification: 65P25, 81C06, 81C07 



I Introduction 

Several years ago Neuberger and Noid |l]j|] presented 
an algorithm and a FORTRAN program for the suc- 
cessive computation of the high-lying eigenvalues and 
eigenfunctions of a time independent Schrodinger or 
Helmholtz equation. They used an inverse power method 
with a Gauss-Seidel procedure for the inversion, and 
solved the problem with finite differences on successively 
finer grids. This program was often used and a two- 
dimensional version thereof @ was adapted for the case 
of Laplace operators with homogeneous boundary condi- 
tions 

HI- The case of Dirichlet conditions works very 
well and requires minimal adjustments. This is not the 
case for Neumann conditions. Computation times in- 
crease by orders of magnitude compared to the Dirichlet 
case 0]; this is not surprising as the treatment of an 
irregular boundary, and particularly of corners, is very 
cumbersome. 

The need for such programs arises both in acoustic || 
and earthquake research J5j, as well as for other wave 
phenomena. For example, if we want to make statistics 
of eigenvalues, say for acoustic systems ||, we need large 



numbers of states and therefore efficient algorithms, fn 
particular, geometries whose high-frequency limit (ray 
dynamics) is chaotic, are of great interest. The start- 
ing point in this new field called acoustic chaos is that 
the time-independent wave equation (Helmholtz equa- 
tion) is the same for different systems such as: quantum 
billiards |s|JTo|] , membranes |ll[] and flat microwave cav- 
ities [|I2|-[14|]. Thus the statistical fluctuation measures 
developed in nuclear physics |l5],|l6| have been applied to 
those systems and to a wide variety of more complicated 
systems such as: Chladni's plates jl7|,[l8), quartz crys- 
tals aluminum blocks pO 21|, quantum dots [ p2[ , 
quantum corrals |2^,Q , waves in a ripple tank |25|] , elas- 
tic media with ray splitting ]2l|, microwave cavities with 
ray splitting |2^] amog others. As rather fine details of 
the boundary of those systems are believed to be impor- 
tant, a good representation of the boundary conditions 
is essential. Similar arguments will hold if we wish to 
study the effect of obstacles inside a cavity or in the old 
Tenochtitlan lake bed || . 

We shall use the finite element method (FEM). It is 
based on the minimization of the functional: 
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T [*] = / (W) 2 ds ~ k 2 [ * 2 ds, 
JR jr. 



(1) 



where i? is a two dimensional region and ip is the wave 
function. 

The possibility of solving mixed boundary conditions 
will be important in systems with mirror symmetries, in 
which we may solve each non-symmetric part using mixed 
boundary conditions. An irregular boundary as well as 
corners can also easily be implemented. The FEM for- 
malism is based on minimizing the functional (|]J) not 
in the complete Hilbert space but only in a sub-space 
spanned by a finite set of piece- wise linear functions, typ- 
ically defined as pyramids over hexagonal cells. We call 
an element the triangles that form this cell. 

The finite difference method becomes very cumber- 
some for boundary conditions involving normal deriva- 
tives at irregular boundaries. This is particularly trou- 
blesome for Neumann conditions for which it leads to 
poor convergence. Novaro et al. ( 0) found in a partic- 
ular case that computations would be two orders of mag- 
nitude slower for Neumann conditions than for Dirichlet 
conditions. 

The minimization of the functional ([!]) on the subspace 
of Hilbert space and in the non-orthogonal basis men- 
tioned above yields the generalized eigenvalue equation 



Ax = XBx 7 



(2) 



with Aij = J R \/ipi ■ \7^jds and Bij = J R i/jitjjjds, where 
ipi and ipj denote the functions defined around the node i 
and j of the hexagonal grid respectively. These functions 
for simplicity, are taken to be linear 



(3) 



with a^, b\ and c\ constants to be determined for each 
triangle of the i-th hexagonal element. The trial func- 
tions are then defined as pyramids of unit height over 
each hexagonal cell, i.e. one function corresponding to 
each node of the grid. The simplicity of a piecewise lin- 
ear basis gives as result a non-orthogonal basis. The Neu- 
mann boundary conditions are obtained by allowing vari- 
ations of the trial functions on the boundary. The Dirich- 
let boundary conditions can be obtained by putting the 
trial functions to zero in the desired part of the bound- 
ary. We refer to the literature for a general and deeper 
discussion of the finite element formalism 28 31 1. 

If the dimension N of the matrices A and B is small 
enough that they can be diagonalised, standard tech- 
niques for non-orthogonal bases may be used. But in 
a typical application, the dimensions are of several thou- 
sand to tenthousands. Yet we are only interested in 
a fairly small number of eigenvalues and eigenfunctions 
near the low end of the spectrum. We thus have to use 
some method that makes explicit use of the sparseness of 
the matrices A and B. We shall see in the next section, 



that a combination of the inverse power and Gauss-Seidel 
methods proposed by Neuberger and Noid can be 
generalised to solve Eq. (||). 

Thus in the following section we show how the inverse 
power method can be applied when a non-orthogonal ba- 
sis is used. Next we discuss how to implement this for 
finite differences as well as a number of tricks that can be 
used to accelerate the numerical procedure and comment 
on the performance of the program. In section IV we 
apply the program to the acoustic stadium, analyze the 
resulting spectra and the corresponding states. We give 
a correction to the spectral density based on an analysis 
of the equations infinite elements for the exactly solvable 
rectangle discussed in the appendix, Finally we present 
some conclusions. 

II The inverse power method in a 
non-orthogonal basis 



We shall transform the Eq. (g) 
with B^ 1 to the form 

B~ x Ax. = Ax 



by left multiplication 



(4) 



As usual a power (B~ 1 A)~ M = (A~ 1 B) M applied to 
an arbitrary initial vector x° will successively select the 
lowest eigenvector corresponding to eigenvalue Ax as it 
will appear with a power (l/Ai) M . 

The problem now seems to be that A -1 is no longer a 
sparse matrix, but this can be averted in the procedure 
of applying A~ 1 B. Thus we need to know the product 



y^A^Bx . 

In order to obtain this product we define x° 
giving for Eq. (||) 

y = A- 1 * 
Multiplying by A by the left we obtain 
Ay = x° 



(5) 
Bx° 

(6) 
(7) 



that can be solved alternatively by Jacobi or Gauss-Seidel 
procedures which will work well as all matrices involved 
continue to be sparse. Summarising: the Gauss-Seidel 
Method can be utilized because we do not need the in- 
verse matrix A~ x but the product A~ x yp . 

Up to here we have only specified how to obtain the 
lowest state; for excited states the usual procedure is to 
orthogonalize the space in which we carry out the varia- 
tion on all states that have already been calculated. Here 
again the non-orthogonality of our basis has to be taken 
into account; indeed, the eigenstates of the Laplace op- 
erator are orthogonal and we have to derive the conse- 
quences this has for eigenvectors in our non-orthogonal 
basis. If and are eigenstates of the Laplace oper- 
ator and are expanded as with coefficients a l m ; I = i,j 
that form vectors x; we can readily check that 
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/ * i B* i ds = 5i j (8) 

JR. 

implies 

x|Bx 3 - % (9) 

and vice-versa. Thus we replace the usual orthogonaliza- 
tion by what we may call a i?-orthogonalization, i.e. we 
require the new vectors to be orthogonal on Bxj (where 
j indicates the calculated eigenstates) thus guaranteeing 
the validity of Eq. (^) and by consequence Eq. (||) . 

Ill Implementation of finite elements and 
convergence 

Once we know the grid and the matrix elements of A 
and B both of which will be evaluated in the end of this 
section, we are in principle ready to write our program. 
As usual a program consists in part of an efficient algo- 
rithm which we have presented, and in part of a bag of 
semi-empirical tricks that tend to repeat themselves in 
different guises again and again. The efficiency of the 
latter is quite problem-dependent and the corresponding 
parameters must be adjusted to optimize operation of the 
program in every case. We shall give some recommenda- 
tions that ought to work reasonably, but we urge the user 
to fiddle around with these parameters. 

When using Neumann boundary conditions, the lowest 
eigenvalue is zero and its corresponding eigenfunction is 
a constant. Yet the Gauss-Seidel inversion requires pos- 
itive definite eigenvalues. This is obtained by adding an 
arbitrary constant C multiplied by B to the Laplace op- 
erator. If we choose this constant large we will need few 
Gauss-Seidel iterations as the operator is near diagonal. 
On the other hand, we will loose precision and conver- 
gence speed in the inverse power process as neighbour- 
ing eigenvalues will have their inverse very close to each 
other. A good balance seems to be to choose the constant 
smaller than, but of the order of, the largest eigenvalue 
which we want to obtain. In the case of Dirichlct condi- 
tions this parameter can also be introduced as a means 
of improving convergence exclusively. Adequate choices 
of this parameter may improve computation time by a 
factor of ~ 2. 

Superconvergence, or over-relaxation, is another stan- 
dard tool to improve convergence. We introduce a fac- 
tor 1 < a < 2 and at each step of any iteration from 
x„ — > x n+ i we replace x n+1 by x„ + a(x n+ i - x„) thus 
correcting a little more than the iteration warrants. For 
the Gauss-Seidel iterations values of a w 1.5 yield im- 
provements in computing time of the order of ~ 2, which 
is consistent with what we found for finite differences. 
But a careful analysis for the stadium shows that at least 
in this case for a quite precise value of a = 1.75 we find 
an acceleration by a factor of ~ 3. On the other hand, 
in the power iteration, improvements are not very signif- 
icant and only values of a near 1 seem acceptable. We 



do not use superconvergence in this context, but again 
we must warn that it might be very significant in other 
cases. All these parameters were carefully explored by 
Mendez 0. 

Another option is the stepwise introduction of finer 
grids, with interpolated results from the rougher grid 
results as starting point. This idea was extensively ex- 
ploited in the finite difference programs of Neuberger and 
Noid jl]]2| and gave excellent results both shortening the 
computation time by giving a good trial function on the 
finest grid and allowing a considerable improvement of 
the eigenvalue upon extrapolation. Unfortunately these 
advantages diminish even in their case as we go to higher 
states, for which only the finest grid is acceptable (finer 
ones would yield too large matrices). For this reason 
we have not implemented these procedures for finite ele- 
ments at this point. 

Now we return with the problem of establishing the 
grid that defines our finite elements and of calculating the 
matrix elements of A and B. For this we use a standard 
method |28| -|nj summarised in the following steps: 

1. We immerse the region R in a quadratic grid. The 
triangles of the elements are defined using the sides 
of the squares and in addition one of the diagonals. 

2. We redefine the grid and triangles along the bound- 
ary by considering all points that lie just outside 
our contour. We then move the exterior points 
along the edges of the squares or the selected di- 
agonal, to the boundary, so as not to change the 
topology of the grid and its connections. 

The corresponding integrals are evaluated using a 
change of variable with a linear transformation. The 
transformation can be found solving a linear system of 
3 equations. Evaluating the constants of Eq. (^) and the 
Jacobian of the transformation (in this case the one half 
area of the triangle, because the transformation is linear) 
we obtain for the integrals: 

A ij =Y J S{^k){a k i a k j +b k i b>;) (10) 

k 

and 

f E fc 2S(A fc )i,^i, 
Bn = { (11) 

I E fc 2S(A fe )^; = j 

Here 5(A&) is the area of triangle A& and the sum is 
made over the number of triangles common to both ele- 
ments i and j . 

The routines consist of two main parts. The first com- 
putes all elements different from zero. The elements are 
put in two matrices of N x 7 to use the sparseness of 
matrices. In order to retain the original positions of each 
element we construct an index matrix. The second part 
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solves the generalized system of equations using the in- 
verse power iteration or a standard diagonalization free- 
ware routine |3!| , if the dimensions of the matrices are 
small (less than « 5000 x 5000). We obtain w 500 eigen- 
values with error less than 5% of the average spacing 
between levels. 

The input for the first routine is a set of points lo- 
calised on the border of the region of interest and a list 
of intervals of this enumeration in which we want Dirich- 
let boundary conditions. We assume Neumann boundary 
conditions in the rest of the border. 

For the second routine the inputs are the number of 
eigenvectors, the tolerance for the Gauss Seidel and in- 
verse power iterations, the constant C, and the supercon- 
vergence factor a. The output consists of two files with 
the corresponding eigenvalues and eigenvectors. 

IV Application to the acoustic stadium 

By way of example we apply the program to some partic- 
ular two-dimensional region R; such as the Bunimovich 
stadium which is completely chaotic in classical mechan- 
ics. This geometry consists of two semicircles of radius r, 
joined by two straight lines with length 21 |4Q] . We will 
apply pure Neumann boundary conditions. On a grid of 
~ 3000 points we obtain, by the inverse power iteration 
method described above, 200 eigenvectors and eigenval- 
ues with a CPU time of approximately 200 sec. on an 
ALPHA Work-Station. 

To analised a discrete sequence of eigenvalues we first 
define the cumulative level density or staircase function 



A' 



N(E)='£®(E-E i ) 



(12) 



i=i 



and its derivative p{E) (the level density), 



p(E)=J2S(E-E t ). 



(13) 



Here O and S are the Heaviside and Dirac delta functions 
respectively. The staircase function is usually divided in 
a smooth part plus a fluctuating part: 



N{E)=N(E)+N fluct (E). 



(14) 



The cumulative level density obtained by the finite el- 
ement method for the stadium with Neumann boundary 
conditions is plotted in Fig. la. For comparison the av- 
erage level density 



N 



Weyl 



(E) = {AE + PVE 



(15) 



obtained from the Weyl formula |^l[ is depicted. Here A 
is the area of the billiard, P is the length of its perime- 
ter and k is a constant that contains information on the 
topological nature of the billiard and the curvature of 
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FIG. 1. a) Cummulative level density for a quarter of Buni- 
movich stadium with Neumann boundary conditions and dis- 
cretised with 2791 elements. The dashed curve correspond 
to the theoretical result given by the Weyl formula and the 
dotted curve is the Weyl formula corrected by the polyno- 
mial from the equations in finite elements for the square, b) 
Cummulative level density for a 50 x 50 square with periodic 
boundary conditions. As reference the theoretical prediction 
for the square is also plotted. 



its boundary. From this figure we can see that the fi- 
nite element method gives a systematic deviation of the 
eigenvalues || . This systematic deviation in the cumula- 
tive level density is due to the discretisation of the finite 
element. To show this we calculate in the appendix the 
eigenvalues A rl , TO for the equations in finite elements for 
an a x a square with periodic boundary conditions. The 
final equation is 



4 — 2 (cos (k x ) + cos (k v )) 



+ g (cos (k x ) + cos (k y ) + cos (k x + k y )) 



(16) 



where k x = — and k v = — are the wave number on the 
x and y directions, a is the size of the side and n,m — 
0,±1,±2,... 

In Fig. lb we show the cumulative level density ob- 
tained from Eq. ( |l6| ) for a 50 x 50 square. We also 
plotted the cumulative level density for the exact square 
(oc k x + ky). The one coming from the diagonalisation is 
indistinguishable from the obtained by the Eq. (|l^) . The 
systematic deviation observed in the square with periodic 
boundary conditions will be used to correct the Weyl for- 
mula for arbirary-shaped billiards and arbitray boundary 
conditions. To do this we calculate the difference between 
the fits for both cumulative level densities -Eq. ([lif)- and 
the exact square. The polynomial for the difference is 
A(4.60055 - 13.75335S - 10.44418£; 2 + 0.45601£ 3 ) /2500, 
where A is the area of the billiard. If this polynomial is 
added to the Weyl prediction, the resulting curve yields 
good agreement up to ~ 400 eigenvalues in the stadium 
(See Fig. la). 
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FIG. 2. Nearest-Neighbour spacing distribution for the 
GOE (solid line), and for the Poisson sequence (dotted line). 
The one obtained for the stadium (200 eigenvalues for each 
symmetry-NN,DD,ND,DN- of a quarter of the Bunimovich 
stadium with Neumann boundary conditions) is depicted by 
the histogram. 

We shall study the fluctuation properties of the spec- 
trum. In order to do this for a typical sequence of eigen- 
values, it is necessary to suppress the secular variation. 
This "unfolding" of the levels is done through the map- 
ping 



Bk .— > E[ = N(Ei 



PWigner(S) = -sexp(--S 



s > 



FIG. 3. The width er(fc) of the fcth-neighbour spacing dis- 
tribution p(k; Sfe) as a function of k for the GOE (solid line), 
the Poisson (dotted line) and the acoustic stadium (squares). 



Z 2 (L) 



(17) 



The effect of such mapping on the original sequence 
is that the new sequence has on the average a constant 
spacing equal to one. Although we should use the Weyl 
formula corrected by the polynomial, we can use a poly- 
nomial fit N(E) for N(E) that takes into account the 
systematic deviation due to the discretisation. We can 
then calculate different statistical measures of the fluctu- 
ating part of the spectrum. The first statistic we shall use 
is the nearest-neighbour spacing distribution p(s) where 
s = E'^-y — E^, which gives information on the short range 
correlations. The p(s) for the stadium is given in Fig. 2. 
Notice that it agrees with the values predicted by the 
gaussian orthogonal ensemble (GOE) typical for chaotic 
systems. For completeness, the Poisson case, typical for 
integrable systems, is also depicted. The spectrum of the 
acoustic Bunimovich stadium shows p(s — > 0) — * 0. This 
behaviour is called level repulsion. In fact the nearest- 
neighbour spacing distribution agrees with 



(18) 



know in spectral statistics as the Wigner surmise, which 
is very close to the GOE prediction. 

We can also define the /cth-neighbour spacing distri- 
bution p(k;sk)- Now Sk = E' i+k _ i ,— E[, so that p(s) = 
p(0; so = s). It is well known W& that these distribu- 
tions tend to a normal distribution as k grows. Since 
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FIG. 4. The number variance Y?(L) as a function of the 
length L of the interval for the same cases as in Fig. 3. 



(sfe) = k+ 1, the only relevant parameter left is the width 
cr(fc) of the distribution. In Fig. 3 we show a(k) for the 
stadium, GOE and Poisson cases. 

The correlation coefficient between adjacent spacings is 
another short range fluctuation measure. For the acous- 
tic stadium we obtained C = —0.26 near to the GOE 
value (Cgoe = —0.27) and far from the Poisson value 

(C Poisson 0). 

Another commonly used statistic is the number vari- 
ance E 2 (L), defined as the second moment of the num- 
ber of levels v{L) within an interval of length L, and 
given for the stadium, GOE and the Poisson cases in 
Fig. 4. For the stadium and GOE cases and for large L, 
S 2 (L) w ln(£) indicating a very rigid sequence. 

The number variance depends exclusively on the two- 
point function. We consider further moments 



V k {L) = ([ V {L)-(v{L))] h 



(19) 
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FIG. 6. The excess 72(1/) -Eq. (21)- as a function of the 
length L of the interval for the same cases as in Fig. 3. 

that depend on higher correlations. To emphasise the 
three- or four-point properties ]32|] , it is useful to consider 
the skewness 

7l (L) = S 3 (i)x [Y?(L)Y V2 (20) 

and the excess 

72 (ri) = £ 4 (L) x [£ 2 (L)f 2 -3. (21) 
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time (level density) 

FIG. 7. The power spectrum \C(t)\ 2 -Eq. (22)- as a func- 
tion of the dimensionless variable t, for the same cases as in 
Fig. 3. We eliminated the divergence at the origin due to the 
finite range of the spectrum. 

depends exclusively on the two-point function ]33[ |. The 
short-range part of Pit) gives specific information con- 
cerning the long-range stiffness. Numerical values are 
given in Fig. 7. 

The eigenfunctions for the rectangle were also success- 
fully compared with the exact solution. In Fig. 8 we 
show two eigenfunctions of the stadium: one of them 
with Dirichlet boundary conditions which shows typical 
feature of the whispering gallery states j$4j , and the last 
one, with Neumann boundary conditions shows a near 
bouncing ball state. There are other kinds of features 
in different eigenfunctions, like scars [ j35|]36| , and others 
that resemble noise |57|,[}8]]. Some which have been re- 
produced by the autors with this algorithm can be seen 
in ref. 0. 

Finally, we performed all the calculation on a quar- 
ter stadium and used Neumann (N) or Dirichlet (D) 
boundary conditions on the symmetry lines. This im- 
plies that we studied the symmetric or antisymmetric 
solutions with respect to both reflection symmetries of 
the stadium. The full solution, shown in the figures, is 
recovered making the corresponding reflections respect 
each symmetry axis. 

V Conclusions 



The numerical values obtained for the stadium are 
shown in Figs. 5 and 6. 

In many instances the Fourier transform of the spectra 
has also proven useful. It contains the same information 
as the spectrum itself. On the other hand the power 
spectrum: 



\c{t)Y 



i 

2^ 



-2i7rEt 



p{E)dE 



(22) 



We have presented an algorithm based on the finite el- 
ement method that computes eigenfunctions and eigen- 
values of the two-dimensional Hclmholtz equation with 
mixed Neumann and Dirichlet boundary conditions. The 
algorithm is divided in two parts: one that computes the 
matrix elements and another that diagonalizes the gen- 
eralised eigenvalue equation using the inverse power and 
Gauss-Seidel Methods. The Gauss-Seidel method runs 
efficiently if we use over-relaxation where the gain in com- 
puting time peaks at a factor of ss 3 around the value 1.75 
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FIG. 8. Eigenfunctions of Bunimovich stadium calculated 
by finite element method for a quadrant and reflecting with 
respect to both axis: a) Dirichlet boundary conditions with 
a relation l/r = 1/2. This figure show a whispering gallery 
state; b) Neumann boundary conditions with l/r = I. Typical 
state scarred by a near "bouncing-ball" orbit. 



Appendix: Eigenvalues for the equations in 
finite elements. 

In this appendix we deduce the equations in finite el- 
ements for a square with periodic boundary conditions. 
If we denote by 



_ i(fc x n+fc„m) 



the amplitude in the grid point (n,m), the Eq. 
the bulk element is given by 



(23) 
I) for 
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-ik x 



— e - lk -» (24) 
12 v ; 



Here k x — — and k v — — . We also assumed that the 
size of the grid is one and that for the bulk A^i = 4, 
Ai,j±i = Ai±\ t j = — 1, Bij±i = Bi±ij = Bi±i.j±i = j2 
and Bi i = The final form for the eigenvalue equation 
in finite elements is given by 

4—2 (cos(k x ) + cos(ky)) 

,\„.,„ ( - + - (cos(k x ) + cos(ky) + cos(k x + k y ))j . (25) 



for the superconvergence (over-relaxation) . A systematic 
error in frequencies was found. This error is due to the 
discretisation and can be estimated by using the eigen- 
values of the equations in finite elements. The programs 
were used to calculate the eigenvalues of the acoustic sta- 
dium. The fluctuation measures of the eigenvalues were 
compared with the random matrix predictions. 

The algorithm is very useful in diverse branches of 
wave physics. The program can obtain eigenfunctions 
and eigenvalues of two-dimensional acoustic cavities, two- 
dimensional microwave cavities, membranes, quantum 
billiards, elastic valleys in certain approximations, etc., 
and can be readily generalised to other problems. 
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